
    GetContrast = function(data, alpha, beta){
        
        nb = nrow(data)
        l0 = data[,c("l0.1","l0.2")]
        a0 = data[,"a0"]
        theta  = exp(l0 %*% t(alpha))
        phi123 = exp(l0 %*% t(beta[1:3,]))
        phi45  = expit(l0 %*% t(beta[4:5,]))
        phi    = cbind(phi123, phi45)
        p.y.cond = t(sapply(1:nb, function(i) GetProb(theta[i,],phi[i,])))
        colnames(p.y.cond) = c("p000","p001","p010","p011","p100",
                               "p101","p110","p111")
        p.l1 = sapply(1:nb, function(i) phi45[i,2 - a0[i]])
        
        EY.00 = mean(p.y.cond[,"p010"] * phi45[,2] + 
                         p.y.cond[,"p000"] * (1-phi45[,2]))
        EY.11 = mean(p.y.cond[,"p111"] * phi45[,1] + 
                         p.y.cond[,"p101"] * (1-phi45[,1]))
        EY.11.00 = EY.11 / EY.00
        
        return(c(EY.00, EY.11, EY.11.00))
        
    }